#ifndef ELEMENT_H
#define ELEMENT_H

#include "Dofs.h"
#include "GaussianPoint.h"
#include <features.h>

/**
 * @brief 有限元计算过程中的基类单元，用于派生出各种不同单元；
 *        Element是TemplateElement的缩写，表示模板单元或者成为referernce element;
 */
template<unsigned int DIM>
class Element
{

public:
    /**
     * @brief Construct a new Element object
     * 
     */
    Element() = default;
    /**
     * @brief Destroy the Element object
     * 
     */
    ~Element() = default;
    /**
     * @brief Set the Dofs List object by std::vector<Dofs<DIM> >
     * 
     */
    virtual void SetDofsList(std::vector<Dofs<DIM> >) = 0 ; 
    /**
     * @brief Get the Dofs List object
     * 
     * @param i 指定序号；
     * @return Dofs<DIM>&　自由度的引用 
     */
    virtual Dofs<DIM>& GetDofsList(int i) = 0 ;
    /**
     * @brief 返回本单元里第ｉ个自由度对应的全局编号；
     * 
     * @param i 指定序号；
     * @return int 全局编号；
     */
    virtual int NdIdx(int i)  = 0 ;
    /**
     * @brief 返回本单元对应第i个高斯点坐标信息；
     * 
     * @param i 指定序号
     * @return double* 存储高斯点的数组
     */
    virtual double* GaussionPoint(int i) = 0 ;
    /**
     * @brief 返回本单元对应第i个高斯点权重信息；
     * 
     * @param i 指定序号
     * @return double 高斯点对应权重的数值；
     */
    virtual double GaussionWeight(int i) = 0 ;
    /**
     * @brief 该类型单元对应的基函数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i　基函数序数
     * @return double 该基函数的值；
     */
    virtual double phi(double xi, double eta, int i)= 0 ;
    /**
     * @brief 该类型单元对应的基函数对xi的偏导数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 该偏导数的值
     */
    virtual double phi_xi(double xi, double eta, int i)  = 0;
    /**
     * @brief 该类型单元对应的基函数对eta的偏导数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 该偏导数的值
     */
    virtual double phi_eta(double xi, double eta, int i)  = 0;
    /**
     * @brief 任意单元到该模板单元仿射变换对应的雅克比矩阵行列式在(xi,eta)处的值；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @return double 值
     */
    virtual double det_Jacobi(double xi, double eta) = 0 ;
    /**
     * @brief 基函数对x的偏导数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 值
     */
    virtual double phi_x(double xi, double eta, int i)  = 0 ;
    /**
     * @brief 基函数对y的偏导数；
     * 
     * @param xi 坐标１
     * @param eta 坐标２
     * @param i 序数
     * @return double 值
     */
    virtual double phi_y(double xi, double eta, int i)  = 0 ;
    /**
     * @brief 计算并返回第ｉ个基函数在xi，eta点处的梯度向量；
     * 
     * @param xi 
     * @param eta 
     * @param i 
     * @return std::vector<double> 梯度向量；
     */
    virtual std::vector<double> gradient(double xi,double eta,int i) = 0;
    /**
     * @brief 标准单元中的点在真实单元的x坐标
     * 
     * @param xi 坐标１
     * @param eta 坐标2
     * @return double 值
     */
    
    virtual double Global_x(double xi, double eta)  = 0 ;
    /**
     * @brief 标准单元中的点在真实单元的y坐标
     * 
     * @param xi 坐标１
     * @param eta 坐标2
     * @return double 值
     */
    virtual double Global_y(double xi, double eta)  = 0;
    /**
     * @brief 该单元自由度个数；
     * 
     * @return int 
     */
    virtual int n_Dofs() = 0;
    /**
     * @brief 高斯点个数；
     * 
     * @return int 
     */
    virtual int n_GaussPnt() = 0;
    /**
     * @brief 计算并返回该单元的面积或体积；
     * 
     * @return double 
     */
    virtual double volume() = 0;
    protected:
      /**
       * xi对x的偏导数.xi_x = y_eta/det_jacobi
       *                   = (y1*phi_eta(1) + y2*phi_eta(2) +y3*phi_eta(3))/det_jacobi;
       *
       * @param xi 坐标1
       * @param eta 坐标2
       *
       * @return double
       */
      virtual double xi_x(double xi, double eta) = 0;
          /**
       * xi对y的偏导数.xi_y = -x_eta/det_jacobi
       *                   = (x1*phi_eta(1) + x2*phi_eta(2) +x3*phi_eta(3))/det_jacobi;
       *
       * @param xi 坐标1
       * @param eta 坐标2
       *
       * @return double
       */
      virtual double xi_y(double xi, double eta) = 0;
          /**
       * eta对x的偏导数.eta_x = -y_xi/det_jacobi
       *                   = (y1*phi_xi(1) + y2*phi_xi(2) +y3*phi_xi(3))/det_jacobi;
       *
       * @param xi 坐标1
       * @param eta 坐标2
       *
       * @return double
       */
      virtual double eta_x(double xi, double eta) = 0;
          /**
       * eta对y的偏导数.eta_y = x_xi/det_jacobi
       *                   = (x1*phi_xi(1) + x2*phi_xi(2) +x3*phi_xi(3))/det_jacobi;
       *
       * @param xi 坐标1
       * @param eta 坐标2
       *
       * @return double
       */
    virtual double eta_y(double xi, double eta) = 0;
    
};



/***************************************************************/
/********************** Q1 ELEMENT *****************************/
/***************************************************************/
class Quadrilateral_1_Element:public Element<2>
{
    typedef std::vector<Dofs<2> > DofsList;
    typedef int DofsNum;
public:
    Quadrilateral_1_Element();
    ~Quadrilateral_1_Element();
    virtual void SetDofsList(DofsList);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual double* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual double GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};               
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i)  ;
    virtual double phi_eta(double xi, double eta, int i)  ;
    virtual double det_Jacobi(double xi, double eta)   ;
    virtual double xi_x(double xi, double eta)  ;
    virtual double xi_y(double xi, double eta)   ;
    virtual double eta_x(double xi, double eta)   ;
    virtual double eta_y(double xi, double eta)   ;
    virtual double phi_x(double xi, double eta, int i)  ;
    virtual double phi_y(double xi, double eta, int i)  ;
    virtual std::vector<double> gradient(double xi,double eta,int i) ;
    virtual double Global_x(double xi, double eta) ;
    virtual double Global_y(double xi, double eta)  ;
    virtual double volume(){return 4 * det_Jacobi(0,0);};
    virtual int n_Dofs(){return 4;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    DofsList _dofslist;
    Quadrilateral_Element_GaussianInfo _GaussInfo;
};

Quadrilateral_1_Element::Quadrilateral_1_Element()
{
    _dofslist.resize(4);
    _GaussInfo.LoadGaussianInfo(3);
};

Quadrilateral_1_Element::~Quadrilateral_1_Element(){};
void Quadrilateral_1_Element::SetDofsList(DofsList _list)
{
    _dofslist = _list;

}
Dofs<2>& Quadrilateral_1_Element::GetDofsList(int i)
{
    return _dofslist[i];
}
double Quadrilateral_1_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (xi-1)*(eta-1)/4.0;
    case 2:
	return -(xi+1)*(eta-1)/4.0;
    case 3:
	return (xi+1)*(eta+1)/4.0;
    case 4:
	return -(xi-1)*(eta+1)/4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::phi_xi(double xi, double eta, int i) 
{
    switch(i)
    {
    case 1:
	return (eta-1)/4.0;
    case 2:
	return -(eta-1)/4.0;
    case 3:
	return (eta+1)/4.0;
    case 4:
	return -(eta+1)/4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::phi_eta(double xi, double eta, int i) 
{
    switch(i)
    {
    case 1:
	return (xi-1)/4.0;
    case 2:
	return -(xi+1)/4.0;
    case 3:
	return (xi+1)/4.0;
    case 4:
	return -(xi-1)/4.0;
    }
    return 0;
}

double Quadrilateral_1_Element::det_Jacobi(double xi, double eta) 
{
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Quadrilateral_1_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::eta_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_1_Element::phi_x(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Quadrilateral_1_Element::phi_y(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}

std::vector<double> Quadrilateral_1_Element::gradient(double xi,double eta,int i)
{
    std::vector<double> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Quadrilateral_1_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Quadrilateral_1_Element::Global_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Quadrilateral_1_Element::NdIdx(int i)
{
    return  _dofslist[i - 1].GlobalIndex();
}

/***************************************************************/
/********************** Q2 ELEMENT *****************************/
/***************************************************************/
class Quadrilateral_2_Element:public Element<2>
{
    typedef std::vector<Dofs<2> > DofsList;
    typedef int DofsNum;
public:
    Quadrilateral_2_Element();
    ~Quadrilateral_2_Element(){};
    virtual void SetDofsList(DofsList);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual double* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual double GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};               
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i)  ;
    virtual double phi_eta(double xi, double eta, int i)  ;
    virtual double det_Jacobi(double xi, double eta)   ;
    virtual double xi_x(double xi, double eta)  ;
    virtual double xi_y(double xi, double eta)   ;
    virtual double eta_x(double xi, double eta)   ;
    virtual double eta_y(double xi, double eta)   ;
    virtual double phi_x(double xi, double eta, int i)  ;
    virtual double phi_y(double xi, double eta, int i)  ;
    virtual std::vector<double> gradient(double xi,double eta,int i) ;
    virtual double volume(){return 4 * det_Jacobi(0,0);};
    virtual double Global_x(double xi, double eta) ;
    virtual double Global_y(double xi, double eta)  ;
    virtual int n_Dofs(){return 9;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    DofsList _dofslist;
    Quadrilateral_Element_GaussianInfo _GaussInfo;
};
Quadrilateral_2_Element::Quadrilateral_2_Element()
{
    _dofslist.resize(9);
    _GaussInfo.LoadGaussianInfo(5);
}
void Quadrilateral_2_Element::SetDofsList(DofsList _list)
{
    _dofslist = _list;
}

Dofs<2>& Quadrilateral_2_Element::GetDofsList(int i)
{
    return _dofslist[i];
}

double Quadrilateral_2_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1/4.0)*xi*eta*(1.0-xi)*(1.0-eta);
    case 2:
	return -(1/4.0)*xi*eta*(1.0+xi)*(1.0-eta);
    case 3:
	return (1/4.0)*xi*eta*(1.0+xi)*(1.0+eta);
    case 4:
    return -(1/4.0)*xi*eta*(1.0-xi)*(1.0+eta);
    case 5:
    return -(1/2.0)*(1.0-eta)*eta*(1.0-xi*xi);
    case 6:
    return (1/2.0)*(1.0+xi)*xi*(1.0-eta*eta);
    case 7:
    return (1/2.0)*(1.0+eta)*eta*(1.0-xi*xi);
    case 8:
    return -(1/2.0)*(1.0-xi)*xi*(1.0-eta*eta);
    case 9:
    return (1.0-xi*xi)*(1.0-eta*eta);
    }
    return 0;
}

double Quadrilateral_2_Element::phi_xi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1/4.0)*(1.0-2.0*xi)*eta*(1.0-eta);
    case 2:
	return -(1/4.0)*(1.0+2.0*xi)*eta*(1.0-eta);
    case 3:
	return (1/4.0)*(1.0+2.0*xi)*eta*(1.0+eta);
    case 4:
    return -(1/4.0)*(1.0-2.0*xi)*eta*(1.0+eta);
    case 5:
    return xi * (1.0-eta) * eta;
    case 6:
    return (1/2.0)*(1.0 + 2.0*xi) * (1.0-eta*eta);
    case 7:
    return -xi * (1.0+eta) * eta;
    case 8:
    return -(1/2.0)*(1.0 - 2.0*xi) * (1.0-eta*eta);
    case 9:
    return -2.0 * xi * (1.0 - eta * eta);
    }
    return 0;
}

double Quadrilateral_2_Element::phi_eta(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1/4.0)*xi*(1.0 - xi)*(1.0-2.0*eta);
    case 2:
	return - (1/4.0)*xi*(1.0 + xi)*(1-2.0*eta);
    case 3:
	return (1/4.0)*xi*(1.0 + xi)*(1.0+2.0*eta);
    case 4:
    return - (1/4.0)*xi*(1.0 - xi)*(1+2.0*eta);
    case 5:
    return - (1/2.0)*(1.0 - xi * xi)*(1.0-2.0*eta);
    case 6:
    return - xi * (1.0+xi) * eta;
    case 7:
    return (1/2.0)*(1.0 - xi * xi)*(1.0+2.0*eta);
    case 8:
    return xi * (1.0-xi) * eta;
    case 9:
    return -2.0 * (1.0 - xi * xi) * eta;
    }
    return 0;
}

double Quadrilateral_2_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 9; i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 9; i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}
double Quadrilateral_2_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 9; i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 9; i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_2_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Quadrilateral_2_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}

std::vector<double> Quadrilateral_2_Element::gradient(double xi,double eta,int i)
{
    std::vector<double> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}

double Quadrilateral_2_Element::det_Jacobi(double xi, double eta) 
{
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 1; i <= 9; i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Quadrilateral_2_Element::Global_x(double xi, double eta) 
{
    /*
    double a = 0;
    for (int i = 1; i <= 4; i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }*/
    double a = 0;
    a += _dofslist[0][0] * 0.25*(1 - xi)*(1 - eta);
    a += _dofslist[1][0] * 0.25*(1 + xi)*(1 - eta);
    a += _dofslist[2][0] * 0.25*(1 + xi)*(1 + eta);
    a += _dofslist[3][0] * 0.25*(1 - xi)*(1 + eta);
    return a;
}

double Quadrilateral_2_Element::Global_y(double xi, double eta)
{
    /*double a = 0;
    for (int i = 1; i <= 4; i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }*/
    double a = 0;
    a += _dofslist[0][1] * 0.25*(1 - xi)*(1 - eta);
    a += _dofslist[1][1] * 0.25*(1 + xi)*(1 - eta);
    a += _dofslist[2][1] * 0.25*(1 + xi)*(1 + eta);
    a += _dofslist[3][1] * 0.25*(1 - xi)*(1 + eta);
    return a;
}

int Quadrilateral_2_Element::NdIdx(int i)
{
    return _dofslist[i - 1].GlobalIndex();
}

/***************************************************************/
/********************** Q3 ELEMENT *****************************/
/***************************************************************/
class Quadrilateral_3_Element:public Element<2>
{
    typedef std::vector<Dofs<2> > DofsList;
    typedef int DofsNum;
public:
    Quadrilateral_3_Element();
    ~Quadrilateral_3_Element();
    virtual void SetDofsList(DofsList);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual double* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual double GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};               
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i)  ;
    virtual double phi_eta(double xi, double eta, int i)  ;
    virtual double det_Jacobi(double xi, double eta)   ;
    virtual double xi_x(double xi, double eta)  ;
    virtual double xi_y(double xi, double eta)   ;
    virtual double eta_x(double xi, double eta)   ;
    virtual double eta_y(double xi, double eta)   ;
    virtual double phi_x(double xi, double eta, int i)  ;
    virtual double phi_y(double xi, double eta, int i)  ;
    virtual std::vector<double> gradient(double xi,double eta,int i) ;
    virtual double Global_x(double xi, double eta) ;
    virtual double Global_y(double xi, double eta)  ;
    virtual double volume(){return 4 * det_Jacobi(0,0);};
    virtual int n_Dofs(){return 16;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    DofsList _dofslist;
    Quadrilateral_Element_GaussianInfo _GaussInfo;
};

Quadrilateral_3_Element::Quadrilateral_3_Element()
{
    _dofslist.resize(16);
    _GaussInfo.LoadGaussianInfo(7);
};

Quadrilateral_3_Element::~Quadrilateral_3_Element(){};
void Quadrilateral_3_Element::SetDofsList(DofsList _list)
{
    _dofslist = _list;

}
Dofs<2>& Quadrilateral_3_Element::GetDofsList(int i)
{
    return _dofslist[i];
}
double Quadrilateral_3_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return 81.0/256 * (xi*xi - 1.0/9) * (xi - 1) * (eta*eta - 1.0/9) * (eta -1);
    case 2:
	return -81.0/256 * (xi + 1) * (xi*xi - 1.0/9) * (eta - 1) * (eta*eta - 1.0/9);
    case 3:
	return 81.0/256 * (xi + 1) * (xi*xi - 1.0/9) * (eta + 1) * (eta*eta - 1.0/9);
    case 4:
	return -81.0/256 * (xi*xi -1.0/9) * (xi - 1) * (eta + 1) * (eta*eta - 1.0/9);
    case 5:
	return -243.0/256 * (xi*xi - 1) * (xi - 1.0/3) * (eta * eta - 1.0/9) * (eta - 1);
    case 6:
	return 243.0/256 * (xi*xi - 1) * (xi + 1.0/3) * (eta*eta - 1.0/9) * (eta - 1);
    case 7:
	return 243.0/256 * (xi + 1) * (xi*xi - 1.0/9) * (eta - 1.0/3) * (eta*eta - 1);
    case 8:
	return -243.0/256 * (xi + 1) * (xi*xi - 1.0/9) * (eta + 1.0/3) * (eta*eta -1);
    case 9:
	return -243.0/256 * (xi*xi - 1) * (xi + 1.0/3) * (eta + 1) * (eta*eta - 1.0/9);
    case 10:
	return 243.0/256 * (xi*xi - 1) * (xi - 1.0/3) * (eta + 1) * (eta*eta - 1.0/9);
    case 11:
	return 243.0/256 * (xi*xi - 1.0/9) * (xi - 1) * (eta + 1.0/3) * (eta*eta - 1);
    case 12:
	return -243.0/256 * (xi*xi - 1.0/9) * (xi - 1) * (eta*eta - 1) * (eta - 1.0/3);
    case 13:
	return 729.0/256 * (xi*xi - 1) * (xi - 1.0/3) * (eta*eta - 1) * (eta - 1.0/3);
    case 14:
	return -729.0/256 * (xi*xi - 1) * (xi + 1.0/3) * (eta*eta - 1) * (eta - 1.0/3);
    case 15:
	return 729.0/256 * (xi*xi - 1) * (xi + 1.0/3) * (eta*eta - 1) * (eta + 1.0/3);
    case 16:
	return -729.0/256 * (xi*xi - 1) * (xi - 1.0/3) * (eta*eta - 1) * (eta + 1.0/3);
    }
    return 0;
}

double Quadrilateral_3_Element::phi_xi(double xi, double eta, int i) 
{
    switch(i)
    {
    case 1:
	return 81.0/256 * (eta*eta -1.0/9) * (eta - 1) * (3*xi*xi - 2*xi - 1.0/9);
    case 2:
	return -81.0/256 * (eta - 1) * (eta*eta - 1.0/9) * (3*xi*xi + 2*xi - 1.0/9);
    case 3:
	return 81.0/256 * (eta + 1) * (eta*eta - 1.0/9) * (3*xi*xi + 2*xi - 1.0/9);
    case 4:
	return -81.0/256 * (eta + 1) * (eta * eta - 1.0/9) * (3*xi*xi - 2*xi - 1.0/9);
    case 5:
	return -243.0/256 *(eta*eta - 1.0/9) * (eta - 1) * (3*xi*xi - 2.0/3*xi - 1);
    case 6:
	return 243.0/256 * (eta*eta - 1.0/9) * (eta - 1) * (3*xi*xi + 2.0/3*xi - 1);
    case 7:
	return 243.0/256 * (eta - 1.0/3) * (eta*eta - 1) * (3*xi*xi + 2*xi - 1.0/9); 
    case 8:
	return -243.0/256 * (eta + 1.0/3) * (eta*eta - 1) * (3*xi*xi + 2*xi - 1.0/9);
    case 9:
	return -243.0/256 * (eta + 1) * (eta*eta - 1.0/9) * (3*xi*xi + 2.0/3*xi - 1);
    case 10:
	return 243.0/256 * (eta + 1) * (eta*eta - 1.0/9) * (3*xi*xi - 2.0/3*xi - 1);
    case 11:
	return 243.0/256 * (eta*eta - 1) * (eta + 1.0/3) * (3*xi*xi - 2*xi - 1.0/9); 
    case 12:
	return -243.0/256 * (eta*eta - 1) * (eta - 1.0/3) * (3*xi*xi - 2*xi - 1.0/9);
    case 13:
	return 729.0/256 * (eta*eta - 1) * (eta - 1.0/3) * (3*xi*xi - 2.0/3*xi - 1);
    case 14:
	return -729.0/256 * (eta*eta - 1) * (eta - 1.0/3) * (3*xi*xi + 2.0/3*xi - 1);
    case 15:
	return 729.0/256 * (eta*eta - 1) * (eta + 1.0/3) * (3*xi*xi + 2.0/3*xi - 1);
    case 16:
	return - 729.0/256 * (eta*eta - 1)* (eta + 1.0/3) * (3*xi*xi - 2.0/3*xi - 1);
    }
    return 0;
}

double Quadrilateral_3_Element::phi_eta(double xi, double eta, int i) 
{
    switch(i)
    {
    case 1:
	return 81.0/256 * (xi*xi - 1.0/9) * (xi - 1) * (3*eta*eta - 2*eta - 1.0/9);
    case 2:
	return -81.0/256 * (xi + 1) * (xi*xi - 1.0/9) * (3*eta*eta - 2*eta - 1.0/9);
    case 3:
	return 81.0/256 * (xi + 1) * (xi*xi - 1.0/9) * (3*eta*eta + 2*eta - 1.0/9);
    case 4:
	return -81.0/256 * (xi*xi - 1.0/9) * (xi - 1) * (3*eta*eta + 2*eta - 1.0/9);
    case 5:
	return -243.0/256 * (xi*xi - 1) * (xi - 1.0/3) * (3*eta*eta - 2*eta - 1.0/9);
    case 6:
	return 243.0/256 * (xi*xi - 1) * (xi + 1.0/3) * (3*eta*eta - 2*eta - 1.0/9); 
    case 7:
	return 243.0/256 * (xi + 1) * (xi*xi - 1.0/9) * (3*eta*eta - 2.0/3*eta - 1);
    case 8:
	return -243.0/256 * (xi + 1) * (xi*xi - 1.0/9) * (3*eta*eta + 2.0/3*eta - 1);
    case 9:
	return -243.0/256 * (xi*xi - 1) * (xi + 1.0/3) * (3*eta*eta + 2*eta - 1.0/9);
    case 10:
	return 243.0/256 * (xi*xi - 1) * (xi - 1.0/3) * (3*eta*eta + 2*eta - 1.0/9);
    case 11:
	return 243.0/256 * (xi*xi - 1.0/9) * (xi - 1) * (3*eta*eta + 2.0/3*eta - 1);
    case 12:
	return -243.0/256 * (xi*xi - 1.0/9) * (xi - 1) * (3*eta*eta - 2.0/3*eta - 1);
    case 13:
	return 729.0/256 * (xi*xi - 1) * (xi - 1.0/3) * (3*eta*eta - 2.0/3*eta - 1);
    case 14:
	return -729.0/256 * (xi*xi - 1) * (xi + 1.0/3) * (3*eta*eta - 2.0/3*eta - 1);
    case 15:
	return 729.0/256 * (xi*xi - 1) * (xi + 1.0/3) * (3*eta*eta + 2.0/3*eta - 1);
    case 16:
	return -729.0/256 * (xi*xi - 1) * (xi - 1.0/3) * (3*eta*eta + 2.0/3*eta - 1);
    }
    return 0;
}

double Quadrilateral_3_Element::det_Jacobi(double xi, double eta) 
{
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 1; i <=n_Dofs() ; i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Quadrilateral_3_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_3_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_3_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_3_Element::eta_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_3_Element::phi_x(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Quadrilateral_3_Element::phi_y(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}

std::vector<double> Quadrilateral_3_Element::gradient(double xi,double eta,int i)
{
    std::vector<double> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Quadrilateral_3_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Quadrilateral_3_Element::Global_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Quadrilateral_3_Element::NdIdx(int i)
{
    return _dofslist[i - 1].GlobalIndex();
}



/***************************************************************/
/********************** Q4 ELEMENT *****************************/
/***************************************************************/
class Quadrilateral_4_Element:public Element<2>
{
    typedef std::vector<Dofs<2> > DofsList;
    typedef int DofsNum;
public:
    Quadrilateral_4_Element();
    ~Quadrilateral_4_Element();
    virtual void SetDofsList(DofsList);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual double* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual double GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};               
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i)  ;
    virtual double phi_eta(double xi, double eta, int i)  ;
    virtual double det_Jacobi(double xi, double eta)   ;
    virtual double xi_x(double xi, double eta)  ;
    virtual double xi_y(double xi, double eta)   ;
    virtual double eta_x(double xi, double eta)   ;
    virtual double eta_y(double xi, double eta)   ;
    virtual double phi_x(double xi, double eta, int i)  ;
    virtual double phi_y(double xi, double eta, int i)  ;
    virtual std::vector<double> gradient(double xi,double eta,int i) ;
    virtual double Global_x(double xi, double eta) ;
    virtual double Global_y(double xi, double eta)  ;
    virtual double volume(){return 4 * det_Jacobi(0,0);};
    virtual int n_Dofs(){return 25;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    DofsList _dofslist;
    Quadrilateral_Element_GaussianInfo _GaussInfo;
};

Quadrilateral_4_Element::Quadrilateral_4_Element()
{
    _dofslist.resize(25);
    _GaussInfo.LoadGaussianInfo(9);
};

Quadrilateral_4_Element::~Quadrilateral_4_Element(){};
void Quadrilateral_4_Element::SetDofsList(DofsList _list)
{
    _dofslist = _list;

}
Dofs<2>& Quadrilateral_4_Element::GetDofsList(int i)
{
    return _dofslist[i];
}
double Quadrilateral_4_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return 4.0/9*xi*(xi*xi - 0.25)*(xi - 1)*eta*(eta-1)*(eta*eta - 0.25);
    case 2:
	return 4.0/9*xi*(xi + 1)*(xi*xi - 0.25)*eta*(eta - 1)*(eta*eta - 0.25);
    case 3:
	return 4.0/9*xi*(xi + 1)*(xi*xi - 0.25)*eta*(eta + 1)*(eta*eta - 0.25);
    case 4:
	return 4.0/9*xi*(xi - 1)*(xi*xi - 0.25)*eta*(eta + 1)*(eta*eta - 0.25);
    case 5:
	return -16.0/9*xi*(xi - 0.5)*(xi*xi - 1)*eta*(eta - 1)*(eta*eta - 0.25);
    case 6:
	return 8.0/3*(xi*xi - 1)*(xi*xi - 0.25)*eta*(eta - 1)*(eta*eta - 0.25);
    case 7:
	return -16.0/9*xi*(xi + 0.5)*(xi*xi - 1)*eta*(eta - 1)*(eta*eta - 0.25);
    case 8:
	return -16.0/9*xi*(xi + 1)*(xi*xi - 0.25)*eta*(eta - 0.5)*(eta*eta - 1);
    case 9:
	return 8.0/3*xi*(xi + 1)*(xi*xi - 0.25)*(eta*eta - 1)*(eta*eta - 0.25);
    case 10:
	return -16.0/9*xi*(xi + 1)*(xi*xi - 0.25)*eta*(eta + 0.5)*(eta*eta - 1);
    case 11:
	return -16.0/9*xi*(xi + 0.5)*(xi*xi - 1)*eta*(eta + 1)*(eta*eta - 0.25);
    case 12:
	return 8.0/3*(xi*xi - 1)*(xi*xi - 0.25)*eta*(eta + 1)*(eta*eta - 0.25);
    case 13:
	return -16.0/9*xi*(xi - 0.5)*(xi*xi - 1)*eta*(eta + 1)*(eta*eta - 0.25);
    case 14:
	return -16.0/9*xi*(xi - 1)*(xi*xi - 0.25)*eta*(eta + 0.5)*(eta*eta - 1);
    case 15:
	return 8.0/3*xi*(xi - 1)*(xi*xi - 0.25)*(eta*eta - 1)*(eta*eta - 0.25);
    case 16:
	return -16.0/9*xi*(xi - 1)*(xi*xi - 0.25)*eta*(eta - 0.5)*(eta*eta - 1);
    case 17:
	return 64.0/9*xi*(xi - 0.5)*(xi*xi - 1)*eta*(eta - 0.5)*(eta*eta - 1);
    case 18:
	return -32.0/3*(xi*xi - 1)*(xi*xi - 0.25)*eta*(eta - 0.5)*(eta*eta - 1);
    case 19:
	return 64.0/9*xi*(xi + 0.5)*(xi*xi - 1)*eta*(eta - 0.5)*(eta*eta - 1);
    case 20:
	return -32.0/3*xi*(xi + 0.5)*(xi*xi - 1)*(eta*eta - 1)*(eta*eta - 0.25);
    case 21:
	return 64.0/9*xi*(xi + 0.5)*(xi*xi - 1)*eta*(eta + 0.5)*(eta*eta - 1);
    case 22:
	return -32.0/3*(xi*xi - 1)*(xi*xi - 0.25)*eta*(eta + 0.5)*(eta*eta - 1);
    case 23:
	return 64.0/9*xi*(xi - 0.5)*(xi*xi - 1)*eta*(eta + 0.5)*(eta*eta - 1);
    case 24:
	return -32.0/3*xi*(xi - 0.5)*(xi*xi - 1)*(eta*eta - 1)*(eta*eta - 0.25);
    case 25:
	return 16*(xi*xi - 1)*(xi*xi - 0.25)*(eta*eta - 1)*(eta*eta - 0.25);
    }
    return 0;
}

double Quadrilateral_4_Element::phi_xi(double xi, double eta, int i) 
{
       switch(i)
    {
    case 1:
	return (eta*(- 4*pow(eta,3) + 4*pow(eta,2) + eta - 1)*(- 16*pow(xi,3) + 12*pow(xi,2) + 2*xi - 1))/36.0;
    case 2:
	return (eta*(- 4*pow(eta,3) + 4*pow(eta,2) + eta - 1)*(- 16*pow(xi,3) - 12*pow(xi,2) + 2*xi + 1))/36.0;
    case 3:
	return  (eta*(- 4*pow(eta,3) - 4*pow(eta,2) + eta + 1)*(- 16*pow(xi,3) - 12*pow(xi,2) + 2*xi + 1))/36.0;
    case 4:
	return (eta*(- 4*pow(eta,3) - 4*pow(eta,2) + eta + 1)*(- 16*pow(xi,3) + 12*pow(xi,2) + 2*xi - 1))/36.0;
    case 5:
	return -(2*eta*(- 4*pow(eta,3) + 4*pow(eta,2) + eta - 1)*(- 8*pow(xi,3) + 3*pow(xi,2) + 4*xi - 1))/9.0;
    case 6:
	return -(eta*xi*(8*pow(xi,2) - 5)*(- 4*pow(eta,3) + 4*pow(eta,2) + eta - 1))/3.0;
    case 7:
	return -(2*eta*(- 4*pow(eta,3) + 4*pow(eta,2) + eta - 1)*(- 8*pow(xi,3) - 3*pow(xi,2) + 4*xi + 1))/9.0;
    case 8:
	return -(2*eta*(- 2*pow(eta,3) + pow(eta,2) + 2*eta - 1)*(- 16*pow(xi,3) - 12*pow(xi,2) + 2*xi + 1))/9.0;
    case 9:
	return -((4*pow(eta,4) - 5*pow(eta,2) + 1)*(- 16*pow(xi,3) - 12*pow(xi,2) + 2*xi + 1))/6.0;
    case 10:
	return -(2*eta*(- 2*pow(eta,3) - pow(eta,2) + 2*eta + 1)*(- 16*pow(xi,3) - 12*pow(xi,2) + 2*xi + 1))/9;
    case 11:
	return -(2*eta*(- 4*pow(eta,3) - 4*pow(eta,2) + eta + 1)*(- 8*pow(xi,3) - 3*pow(xi,2) + 4*xi + 1))/9.0;
    case 12:
	return -(eta*xi*(8*pow(xi,2) - 5)*(- 4*pow(eta,3) - 4*pow(eta,2) + eta + 1))/3.0;
    case 13:
	return -(2*eta*(- 4*pow(eta,3) - 4*pow(eta,2) + eta + 1)*(- 8*pow(xi,3) + 3*pow(xi,2) + 4*xi - 1))/9.0;
    case 14:
	return -(2*eta*(- 2*pow(eta,3) - pow(eta,2) + 2*eta + 1)*(- 16*pow(xi,3) + 12*pow(xi,2) + 2*xi - 1))/9.0;
    case 15:
	return -((4*pow(eta,4) - 5*pow(eta,2) + 1)*(- 16*pow(xi,3) + 12*pow(xi,2) + 2*xi - 1))/6.0;
    case 16:
	return -(2*eta*(- 2*pow(eta,3) + pow(eta,2) + 2*eta - 1)*(- 16*pow(xi,3) + 12*pow(xi,2) + 2*xi - 1))/9.0;
    case 17:
	return (16*eta*(- 2*pow(eta,3) + pow(eta,2) + 2*eta - 1)*(- 8*pow(xi,3) + 3*pow(xi,2) + 4*xi - 1))/9.0;
    case 18:
	return (8*eta*xi*(8*pow(xi,2) - 5)*(- 2*pow(eta,3) + pow(eta,2) + 2*eta - 1))/3.0;
    case 19:
	return (16*eta*(- 2*pow(eta,3) + pow(eta,2) + 2*eta - 1)*(- 8*pow(xi,3) - 3*pow(xi,2) + 4*xi + 1))/9.0;
    case 20:
	return (4*(4*pow(eta,4) - 5*pow(eta,2) + 1)*(- 8*pow(xi,3) - 3*pow(xi,2) + 4*xi + 1))/3.0;
    case 21:
	return (16*eta*(- 2*pow(eta,3) - pow(eta,2) + 2*eta + 1)*(- 8*pow(xi,3) - 3*pow(xi,2) + 4*xi + 1))/9.0;
    case 22: 
	return (8*eta*xi*(8*pow(xi,2) - 5)*(- 2*pow(eta,3) - pow(eta,2) + 2*eta + 1))/3.0;
    case 23:
	return (16*eta*(- 2*pow(eta,3) - pow(eta,2) + 2*eta + 1)*(- 8*pow(xi,3) + 3*pow(xi,2) + 4*xi - 1))/9.0;
    case 24:
	return (4*(4*pow(eta,4) - 5*pow(eta,2) + 1)*(- 8*pow(xi,3) + 3*pow(xi,2) + 4*xi - 1))/3.0;
    case 25:
	return 2*xi*(8*pow(xi,2) - 5)*(4*pow(eta,4) - 5*pow(eta,2) + 1);
    }
    return 0;
}

double Quadrilateral_4_Element::phi_eta(double xi, double eta, int i) 
{
        switch(i)
    {
    case 1:
	return (xi*(- 4*pow(xi,3) + 4*pow(xi,2) + xi - 1)*(- 16*pow(eta,3) + 12*pow(eta,2) + 2*eta - 1))/36.0;
    case 2:
	return (xi*(- 4*pow(xi,3) - 4*pow(xi,2) + xi + 1)*(- 16*pow(eta,3) + 12*pow(eta,2) + 2*eta - 1))/36.0;
    case 3:
	return  (xi*(- 4*pow(xi,3) - 4*pow(xi,2) + xi + 1)*(- 16*pow(eta,3) - 12*pow(eta,2) + 2*eta + 1))/36.0;
    case 4:
	return (xi*(- 4*pow(xi,3) + 4*pow(xi,2) + xi - 1)*(- 16*pow(eta,3) - 12*pow(eta,2) + 2*eta + 1))/36.0;
    case 5:
	return -(2*xi*(- 2*pow(xi,3) + pow(xi,2) + 2*xi - 1)*(- 16*pow(eta,3) + 12*pow(eta,2) + 2*eta - 1))/9.0;
    case 6:
	return -((4*pow(xi,4) - 5*xi*xi + 1)*(- 16*pow(eta,3) + 12*eta*eta + 2*eta - 1))/6;
    case 7:
	return -(2*xi*(- 16*pow(eta,3) + 12*pow(eta,2) + 2*eta - 1)*(- 2*pow(xi,3) - pow(xi,2) + 2*xi + 1))/9.0;
    case 8:
	return -(2*xi*(- 4*pow(xi,3) - 4*pow(xi,2) + xi + 1)*(- 8*pow(eta,3) + 3*pow(eta,2) + 4*eta - 1))/9.0;
    case 9:
	return -(eta*xi*(8*pow(eta,2) - 5)*(- 4*pow(xi,3) - 4*pow(xi,2) + xi + 1))/3.0;
    case 10:
	return -(2*xi*(- 4*pow(xi,3) - 4*pow(xi,2) + xi + 1)*(- 8*pow(eta,3) - 3*pow(eta,2) + 4*eta + 1))/9;
    case 11:
	return -(2*xi*(- 16*pow(eta,3) - 12*pow(eta,2) + 2*eta + 1)*(- 2*pow(xi,3) - pow(xi,2) + 2*xi + 1))/9.0;
    case 12:
	return -((4*pow(xi,4) - 5*pow(xi,2) + 1)*(- 16*pow(eta,3) - 12*pow(eta,2) + 2*eta + 1))/6.0;
    case 13:
	return -(2*xi*(- 2*pow(xi,3) + pow(xi,2) + 2*xi - 1)*(- 16*pow(eta,3) - 12*pow(eta,2) + 2*eta + 1))/9.0;
    case 14:
	return -(2*xi*(- 4*pow(xi,3) + 4*pow(xi,2) + xi - 1)*(- 8*pow(eta,3) - 3*pow(eta,2) + 4*eta + 1))/9.0;
    case 15:
	return -(eta*xi*(8*pow(eta,2) - 5)*(- 4*pow(xi,3) + 4*pow(xi,2) + xi - 1))/3.0;
    case 16:
	return -(2*xi*(- 4*pow(xi,3) + 4*pow(xi,2) + xi - 1)*(- 8*pow(eta,3) + 3*pow(eta,2) + 4*eta - 1))/9.0;
    case 17:
	return (16*xi*(- 2*pow(xi,3) + pow(xi,2) + 2*xi - 1)*(- 8*pow(eta,3) + 3*pow(eta,2) + 4*eta - 1))/9.0;
    case 18:
	return (4*(4*pow(xi,4) - 5*pow(xi,2) + 1)*(- 8*pow(eta,3) + 3*pow(eta,2) + 4*eta - 1))/3.0;
    case 19:
	return (16*xi*(- 8*pow(eta,3) + 3*pow(eta,2) + 4*eta - 1)*(- 2*pow(xi,3) - pow(xi,2) + 2*xi + 1))/9.0;
    case 20:
	return (8*eta*xi*(8*pow(eta,2) - 5)*(- 2*pow(xi,3) - pow(xi,2) + 2*xi + 1))/3.0;
    case 21:
	return (16*xi*(- 8*pow(eta,3) - 3*pow(eta,2) + 4*eta + 1)*(- 2*pow(xi,3) - pow(xi,2) + 2*xi + 1))/9.0;
    case 22:
	return (4*(4*pow(xi,4) - 5*pow(xi,2) + 1)*(- 8*pow(eta,3) - 3*pow(eta,2) + 4*eta + 1))/3.0;
    case 23:
	return (16*xi*(- 2*pow(xi,3) + pow(xi,2) + 2*xi - 1)*(- 8*pow(eta,3) - 3*pow(eta,2) + 4*eta + 1))/9.0;
    case 24:
	return (8*eta*xi*(8*pow(eta,2) - 5)*(- 2*pow(xi,3) + pow(xi,2) + 2*xi - 1))/3.0;
    case 25:
	return 2*eta*(8*pow(eta,2) - 5)*(4*pow(xi,4) - 5*pow(xi,2) + 1);
    }
    return 0;
}

double Quadrilateral_4_Element::det_Jacobi(double xi, double eta) 
{
    double J11 = 0;
    double J12 = 0;
    double J21 = 0;
    double J22 = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Quadrilateral_4_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_4_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_4_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Quadrilateral_4_Element::eta_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Quadrilateral_4_Element::phi_x(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Quadrilateral_4_Element::phi_y(double xi, double eta, int i) 
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}

std::vector<double> Quadrilateral_4_Element::gradient(double xi,double eta,int i)
{
    std::vector<double> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Quadrilateral_4_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Quadrilateral_4_Element::Global_y(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Quadrilateral_4_Element::NdIdx(int i)
{
    return _dofslist[i - 1].GlobalIndex();
}



/***************************************************************/
/********************** P1 ELEMENT *****************************/
/***************************************************************/

class Triangular_1_Element:public Element<2>
{
public:
    Triangular_1_Element();
    ~Triangular_1_Element();
    virtual void SetDofsList(std::vector<Dofs<2> >);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual double* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual double GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};    
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i);
    virtual double phi_eta(double xi, double eta, int i);
    virtual double det_Jacobi(double xi, double eta);
    virtual double xi_x(double xi, double eta);
    virtual double xi_y(double xi, double eta);
    virtual double eta_x(double xi, double eta);
    virtual double eta_y(double xi, double eta);
    virtual double phi_x(double xi, double eta, int i);
    virtual double phi_y(double xi, double eta, int i);
    virtual std::vector<double> gradient(double xi,double eta,int i);
    virtual double volume(){return det_Jacobi(0,0) / 2.0;};
    virtual double Global_x(double xi, double eta);
    virtual double Global_y(double xi, double eta);
    virtual int n_Dofs(){return 3;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    std::vector<Dofs<2> > _dofslist;
    Triangular_Element_GaussianInfo _GaussInfo;
};

Triangular_1_Element::Triangular_1_Element()
{
    _dofslist.resize(3);
    _GaussInfo.LoadGaussianInfo(6);
}

Triangular_1_Element::~Triangular_1_Element()
{
};

void Triangular_1_Element::SetDofsList(std::vector<Dofs<2> > _list)
{
    _dofslist = _list;
}

Dofs<2>& Triangular_1_Element::GetDofsList(int i)
{
    return _dofslist[i];
}

double Triangular_1_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1 - xi - eta);
    case 2:
	return xi;
    case 3:
	return eta;
    }
    return 0;
}

double Triangular_1_Element::phi_xi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return -1;
    case 2:
	return 1;
    case 3:
	return 0;
    }
    return 0;
}

double Triangular_1_Element::phi_eta(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return -1;
    case 2:
	return 0;
    case 3:
	return 1;
    }
    return 0;
}

double Triangular_1_Element::det_Jacobi(double xi, double eta)
{
    return (_dofslist[1][0] - _dofslist[0][0])*(_dofslist[2][1] - _dofslist[0][1])
	-(_dofslist[1][1] - _dofslist[0][1])*(_dofslist[2][0]-_dofslist[0][0]);
}

double Triangular_1_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_1_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Triangular_1_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}
std::vector<double> Triangular_1_Element::gradient(double xi,double eta,int i)
{
    std::vector<double> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Triangular_1_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_1_Element::Global_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= 3; i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Triangular_1_Element::NdIdx(int i)
{
    switch(i)
    {
    case 1:
	return _dofslist[0].GlobalIndex();
    case 2:
	return _dofslist[1].GlobalIndex();
    case 3:
	return _dofslist[2].GlobalIndex();
    }
    return 0;
}


/***************************************************************/
/********************** P2 ELEMENT *****************************/
/***************************************************************/

class Triangular_2_Element:public Element<2>
{
public:
    Triangular_2_Element();
    ~Triangular_2_Element();
    virtual void SetDofsList(std::vector<Dofs<2> >);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual double* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual double GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};                      
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i)  ;
    virtual double phi_eta(double xi, double eta, int i)  ;
    virtual double det_Jacobi(double xi, double eta)   ;
    virtual double xi_x(double xi, double eta)  ;
    virtual double xi_y(double xi, double eta)   ;
    virtual double eta_x(double xi, double eta)   ;
    virtual double eta_y(double xi, double eta)   ;
    virtual double phi_x(double xi, double eta, int i)  ;
    virtual double phi_y(double xi, double eta, int i)  ;
    virtual std::vector<double> gradient(double xi,double eta,int i) ;
    virtual double volume(){return det_Jacobi(0,0)/2.0;};
    virtual double Global_x(double xi, double eta) ;
    virtual double Global_y(double xi, double eta)  ;
    virtual int n_Dofs(){return 6;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    std::vector<Dofs<2> > _dofslist;
    Triangular_Element_GaussianInfo _GaussInfo;
};

Triangular_2_Element::Triangular_2_Element()
{       
    _dofslist.resize(6);
    _GaussInfo.LoadGaussianInfo(8);
}

Triangular_2_Element::~Triangular_2_Element() 
{
};

void Triangular_2_Element::SetDofsList(std::vector<Dofs<2> > _list)
{
    _dofslist = _list;
}

Dofs<2>& Triangular_2_Element::GetDofsList(int i)
{
    return _dofslist[i];
}

double Triangular_2_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return (1 - xi - eta)*(2 *(1 - xi -eta) - 1);
    case 2:
	return xi * (2 * xi - 1);
    case 3:
	return eta * ( 2 * eta - 1);
    case 4:
    return 4 * xi * (1 - xi - eta);
    case 5:
    return 4 * xi * eta;
    case 6:
    return 4 * eta *(1 - xi - eta);
    }
    return 0;
}

double Triangular_2_Element::phi_xi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return 1 - 4*(1- xi - eta);
    case 2:
	return 4 * xi - 1;
    case 3:
	return 0;
    case 4:
    return 4 - 8 * xi - 4 * eta;
    case 5:
    return 4 * eta;
    case 6:
    return - 4 * eta;
    }
    return 0;
}

double Triangular_2_Element::phi_eta(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return 1 - 4*(1- xi - eta);
    case 2:
	return 0;
    case 3:
	return 4 * eta - 1;
    case 4:
    return - 4 * xi;
    case 5:
    return  4 * xi;
    case 6:
    return 4 - 4 * xi - 8 * eta;
    }
    return 0;
}

double Triangular_2_Element::det_Jacobi(double xi, double eta)
{
    double J11 = 0,J12 = 0,J21 = 0,J22 = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Triangular_2_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_2_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Triangular_2_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}
std::vector<double> Triangular_2_Element::gradient(double xi,double eta,int i)
{
    std::vector<double> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Triangular_2_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_2_Element::Global_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Triangular_2_Element::NdIdx(int i)
{
    
    return _dofslist[i - 1].GlobalIndex();
}

class Triangular_3_Element:public Element<2>
{
public:
    Triangular_3_Element();
    ~Triangular_3_Element();
    virtual void SetDofsList(std::vector<Dofs<2> >);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual double* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual double GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};                      
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i)  ;
    virtual double phi_eta(double xi, double eta, int i)  ;
    virtual double det_Jacobi(double xi, double eta)   ;
    virtual double xi_x(double xi, double eta)  ;
    virtual double xi_y(double xi, double eta)   ;
    virtual double eta_x(double xi, double eta)   ;
    virtual double eta_y(double xi, double eta)   ;
    virtual double phi_x(double xi, double eta, int i)  ;
    virtual double phi_y(double xi, double eta, int i)  ;
    virtual std::vector<double> gradient(double xi,double eta,int i) ;
    virtual double volume(){return det_Jacobi(0,0)/2.0;};
    virtual double Global_x(double xi, double eta) ;
    virtual double Global_y(double xi, double eta)  ;
    virtual int n_Dofs(){return 10;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    std::vector<Dofs<2> > _dofslist;
    Triangular_Element_GaussianInfo _GaussInfo;
};

Triangular_3_Element::Triangular_3_Element()
{       
    _dofslist.resize(10);
    _GaussInfo.LoadGaussianInfo(6);
}

Triangular_3_Element::~Triangular_3_Element() 
{
};

void Triangular_3_Element::SetDofsList(std::vector<Dofs<2> > _list)
{
    _dofslist = _list;
}

Dofs<2>& Triangular_3_Element::GetDofsList(int i)
{
    return _dofslist[i];
}

double Triangular_3_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return 4.5 * (1 - xi - eta) * (2.0/3 - xi - eta) * (1.0/3 - xi - eta);
    case 2:
	return 4.5 * xi * (xi - 1.0/3) * (xi - 2.0/3);
    case 3:
	return 4.5 * eta * (eta - 1.0/3) * (eta - 2.0/3);
    case 4:
	return 13.5 * xi * (1 - xi - eta) * (2.0/3 - xi - eta);
    case 5:
	return 13.5 * xi * (1 - xi - eta) * (xi - 1.0/3);
    case 6:
	return 13.5 * xi * eta * (xi - 1.0/3);
    case 7:
	return 13.5 * xi * eta * (eta - 1.0/3);
    case 8:
	return 13.5 * (1 - xi - eta) * eta * (eta - 1.0/3);
    case 9:
	return 13.5 * (1 - xi - eta) * (2.0/3 - xi - eta) * eta;
    case 10:
	return 27 * xi * eta * (1- xi - eta);
    }
    return 0;
}

double Triangular_3_Element::phi_xi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return -4.5 * (2.0/3 - xi - eta) * (1.0/3 - xi - eta) - 4.5 * (1- xi - eta) * (1.0/3 - xi - eta) - 4.5 * (1- xi - eta) * (2.0/3 - xi - eta);  
    case 2:
	return 4.5 * (xi - 1.0/3)*(xi -2.0/3) + 4.5 * xi *(xi - 2.0/3) + 4.5 * xi * (xi - 1.0/3);
    case 3:
	return 0;
    case 4:
	return 13.5 * (xi + eta - 1) * (xi + eta - 2.0/3) + 13.5 * xi * (xi + eta - 2.0/3) + 13.5* xi * (xi + eta - 1);
    case 5:
	return 13.5 * (1 - xi - eta) * (xi - 1.0/3) - 13.5 * xi * (xi - 1.0/3) + 13.5 * xi * (1- xi - eta);  
    case 6:
	return 13.5 * eta * (xi - 1.0/3) + 13.5 * xi * eta;
    case 7:
	return 13.5 * eta * (eta - 1.0/3);
    case 8:
	return -13.5 * eta * (eta - 1.0/3);
    case 9:
	return - 13.5 * (2.0/3 - xi - eta) * eta - 13.5 * (1 - xi - eta)* eta;
    case 10:
	return 27 * eta * (1 - xi - eta) - 27 * xi * eta;
    }
    return 0;
}

double Triangular_3_Element::phi_eta(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return -4.5 * (2.0/3 - xi - eta) * (1.0/3 -  xi - eta)  - 4.5 * (1 - xi -eta) * (1.0/3 - xi - eta) - 4.5 * (1 - xi -eta) * (2.0/3 - xi -eta);
    case 2:
	return  0;
    case 3:
	return 4.5 * (eta - 1.0/3) * (eta - 2.0/3) + 4.5 * eta * (eta - 2.0/3) + 4.5 * eta * (eta- 1.0/3);
    case 4:
	return 13.5 * xi * (xi + eta - 2.0/3) + 13.5 * xi * (xi + eta - 1);
    case 5:
	return  -13.5 * xi * (xi - 1.0/3);
    case 6:
	return 13.5 * xi * (xi -1.0/3);
    case 7:
	return 13.5 * xi * (eta - 1.0/3) + 13.5 * xi * eta;
    case 8:
	return -13.5 * eta * (eta - 1.0/3) + 13.5 * (1 - xi -eta) * (eta - 1.0/3) + 13.5 * (1 - xi - eta) * eta;
    case 9:
	return -13.5 * (2.0/3 - xi - eta) * eta - 13.5 * (1 - xi - eta)* eta + 13.5 * (1 - xi - eta) * (2.0/3 - xi - eta);
    case 10:
	return 27 * xi * (1 - xi - eta) - 27 * xi * eta;
    }
    return 0;
}

double Triangular_3_Element::det_Jacobi(double xi, double eta)
{
    double J11 = 0,J12 = 0,J21 = 0,J22 = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Triangular_3_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_3_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_3_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_3_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_3_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Triangular_3_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}
std::vector<double> Triangular_3_Element::gradient(double xi,double eta,int i)
{
    std::vector<double> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Triangular_3_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_3_Element::Global_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Triangular_3_Element::NdIdx(int i)
{
    
    return _dofslist[i - 1].GlobalIndex();
}



class Triangular_4_Element:public Element<2>
{
public:
    Triangular_4_Element();
    ~Triangular_4_Element();
    virtual void SetDofsList(std::vector<Dofs<2> >);
    virtual Dofs<2>& GetDofsList(int i); 
    virtual int NdIdx(int i) ; 
    virtual double* GaussionPoint(int i){return _GaussInfo.GaussianPointInfo(i);};
    virtual double GaussionWeight(int i){return _GaussInfo.GaussianWeightInfo(i);};                      
    virtual double phi(double xi, double eta, int i);
    virtual double phi_xi(double xi, double eta, int i)  ;
    virtual double phi_eta(double xi, double eta, int i)  ;
    virtual double det_Jacobi(double xi, double eta)   ;
    virtual double xi_x(double xi, double eta)  ;
    virtual double xi_y(double xi, double eta)   ;
    virtual double eta_x(double xi, double eta)   ;
    virtual double eta_y(double xi, double eta)   ;
    virtual double phi_x(double xi, double eta, int i)  ;
    virtual double phi_y(double xi, double eta, int i)  ;
    virtual std::vector<double> gradient(double xi,double eta,int i) ;
    virtual double volume(){return det_Jacobi(0,0)/2.0;};
    virtual double Global_x(double xi, double eta) ;
    virtual double Global_y(double xi, double eta)  ;
    virtual int n_Dofs(){return 15;};
    virtual int n_GaussPnt(){return _GaussInfo.GaussianPointNumInfo();};
private:
    std::vector<Dofs<2> > _dofslist;
    Triangular_Element_GaussianInfo _GaussInfo;
};

Triangular_4_Element::Triangular_4_Element()
{       
    _dofslist.resize(15);
    _GaussInfo.LoadGaussianInfo(8);//查表确定用哪个精度
}

Triangular_4_Element::~Triangular_4_Element() 
{
};

void Triangular_4_Element::SetDofsList(std::vector<Dofs<2> > _list)
{
    _dofslist = _list;
}

Dofs<2>& Triangular_4_Element::GetDofsList(int i)
{
    return _dofslist[i];
}

double Triangular_4_Element::phi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return 32.0/3 * (0.25 - xi - eta) * (0.5 - xi - eta) * (0.75 - xi - eta) * (1 - xi - eta);
    case 2:
	return 32.0/3 * xi * (xi - 0.25) * (xi - 0.5) * (xi- 0.75);
    case 3:
	return 32.0/3 * eta * (eta - 0.25) * (eta - 0.5) * (eta  - 0.75);
    case 4:
	return 128.0/3 * xi * (0.5 - xi - eta) * (0.75 - xi - eta) * (1 - xi - eta);
    case 5:
	return 64 * xi * (xi - 0.25) * (0.75 - xi - eta) * (1 - xi - eta);
    case 6:
	return 128.0/3 * xi * (xi - 0.25) * (xi - 0.5) * (1 - xi - eta);
    case 7:
	return 128.0/3 * xi * (xi - 0.25) * (xi - 0.5) * eta;
    case 8:
	return 64 * xi * (xi - 0.25) * eta * (eta - 0.25);
    case 9:
	return 128.0/3 * xi * eta * (eta - 0.25) * (eta - 0.5);
    case 10:
	return 128.0/3 * eta * (eta - 0.25) * (eta - 0.5) * (1- xi - eta);
    case 11:
	return 64 * eta * (eta - 0.25) * (0.75 - xi - eta) * (1 - xi - eta);
    case 12:
	return 128.0/3 * eta * (0.5 - xi - eta) * (0.75 - xi - eta) * (1 - xi - eta);
    case 13:
	return 128 * xi * eta * (0.75 - xi - eta) * (1 - xi - eta);
    case 14:
	return 128 * xi * (xi - 0.25) * eta * (1 - xi - eta);
    case 15:
	return 128 * xi * eta * (eta - 0.25) * (1 - xi - eta);
    }
    return 0;
}

double Triangular_4_Element::phi_xi(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return -32.0/3 * (0.5 - xi - eta) * (0.75 - xi -eta) * (1 - xi - eta) - 32.0/3 * (0.25 - xi - eta) * (0.75 - xi - eta) * (1 - xi - eta) - 32.0/3 * (0.25 - xi - eta) * (0.5 - xi - eta) * (1 - xi - eta) - 32.0/3 * (0.25 - xi - eta)*(0.5 - xi - eta) * (0.75 - xi - eta); 
    case 2:
	return 32.0/3 * (xi - 0.25) * (xi - 0.5) * (xi - 0.75) + 32.0/3 * xi * (xi - 0.5) * (xi - 0.75) + 32.0/3 * xi * (xi - 0.25) * (xi - 0.75) + 32.0/3 * xi * (xi - 0.25) * (xi - 0.5);
    case 3:
	return 0;
    case 4:
	return 128.0/3 * (0.5 - xi - eta) * (0.75 - xi - eta) * (1 - xi - eta) - 128.0/3 * xi * (0.75 - xi - eta) * (1 - xi -eta) - 128.0/3 * xi * (0.5 - xi - eta) * (1 - xi - eta) - 128.0/3 * xi * (0.5 - xi - eta)* (0.75 - xi - eta);
    case 5:
	return 64 * (xi - 0.25) * (0.75 - xi - eta) * (1 - xi - eta) + 64 * xi * (0.75 - xi - eta) * (1 - xi - eta) - 64 * xi * (xi - 0.25) * (1 - xi - eta) - 64 * xi * (xi - 0.25) * (0.75 - xi - eta);
    case 6:
	return 128.0/3 * (xi - 0.25) * (xi - 0.5) * (1 - xi - eta) + 128.0/3 * xi * (xi - 0.5) * (1 - xi - eta) + 128.0 /3 * xi * (xi - 0.25) * (1 - xi - eta) - 128.0/3 * xi * (xi - 0.25) * (xi - 0.5); 
    case 7:
	return 128.0/3 * (xi - 0.25) * (xi - 0.5) * eta + 128.0/3 * xi * (xi - 0.5) * eta + 128.0/3 * xi * (xi - 0.25) * eta;
    case 8:
	return 64 * (xi - 0.25) * eta * (eta - 0.25) + 64 * xi * eta * (eta - 0.25);
    case 9:
	return 128.0/3 * eta * (eta - 0.25) * (eta - 0.5);
    case 10:
	return -128.0/3 * eta * (eta - 0.25) * (eta - 0.5);
    case 11:
	return -64 * eta * (eta - 0.25) * (1 - xi - eta) - 64 * eta * (eta - 0.25) * (0.75 - xi - eta);
    case 12:
	return -128.0/3 * eta * (0.75 - xi - eta) * (1 - xi - eta) - 128.0/3 * eta * (0.5 - xi - eta) * (1 - xi - eta) - 128.0/3 * eta * (0.5 - xi - eta) * (0.75 - xi - eta);
    case 13:
	return 128 * eta * (0.75 - xi - eta) * (1 - xi - eta) - 128 * xi * eta * (1 - xi - eta) - 128 * xi * eta * (0.75 - xi - eta);
    case 14:
	return 128 * eta * (xi -0.25) * (1 - xi - eta) + 128 * xi * eta * (1 - xi - eta) - 128 * xi * eta * (xi - 0.25);
    case 15:
	return 128 * eta * (eta - 0.25) * (1 - xi - eta) -  128 * xi * eta * (eta - 0.25);
    }
    return 0;
}

double Triangular_4_Element::phi_eta(double xi, double eta, int i)
{
    switch(i)
    {
    case 1:
	return -32.0/3 * (0.5 - xi - eta) * (0.75 - xi - eta) * (1 - xi - eta) - 32.0/3 * (0.25 - xi - eta) * (0.75 - xi - eta) * (1 - xi - eta) -  32.0/3 * (0.25 - xi - eta) * (0.5 - xi - eta) * (1 - xi - eta) - 32.0/3 * (0.25 - xi - eta) * (0.5 - xi - eta) * (0.75 - xi - eta); 
    case 2:
	return  0;
    case 3:
	return 32.0/3 * (eta - 0.25) * (eta - 0.5) * (eta - 0.75) + 32.0/3 * eta * (eta - 0.5) * (eta - 0.75) + 32.0/3 * eta * (eta - 0.25) * (eta - 0.75) + 32.0/3 * eta * (eta - 0.25) * (eta - 0.5);
    case 4:
	return - 128.0/3 * xi * (0.75 - xi - eta) * (1 - xi - eta) - 128.0/3 * xi * (0.5 - xi - eta) * (1 - xi - eta) - 128.0/3 * xi * (0.5 - xi - eta) * (0.75 - xi - eta);
    case 5:
	return  -64 * xi * (xi - 0.25) * (1 - xi - eta) -  64 * xi * (xi - 0.25) * (0.75 - xi - eta);
    case 6:
	return -128.0/3 * xi * (xi - 0.25) * (xi - 0.5);
    case 7:
	return 128.0/3 * xi * (xi - 0.25) * (xi - 0.5);
    case 8:
	return 64 * xi * (xi - 0.25) * (eta - 0.25) + 64 * xi * (xi - 0.25) * eta;
    case 9:
	return 128.0/3 * xi * (eta - 0.25) * (eta - 0.5) + 128.0/3 * xi * eta * (eta - 0.5) + 128.0/3 * xi * eta * (eta - 0.25);
    case 10:
	return 128.0/3 * (eta - 0.25) * (eta - 0.5) * (1 - xi - eta) + 128.0/3 * eta * (eta - 0.5) * (1 - xi - eta) + 128.0/3 * eta * (eta - 0.25) * (1 - xi - eta) - 128.0/3 * eta * (eta - 0.25) * (eta - 0.5);
    case 11:
	return 64 * (eta - 0.25) * (0.75 - xi - eta) * (1 - xi - eta) + 64 * eta * (0.75 - xi - eta) * (1 - xi - eta) -64 * eta * (eta - 0.25) * (1 - xi - eta) - 64 * eta * (eta - 0.25) * (0.75 - xi - eta);
    case 12:
	return 128.0/3 * (0.5 - xi - eta) * (0.75 - xi - eta) * (1 - xi - eta) - 128.0/3 * eta * (0.75 - xi - eta) * (1- xi - eta) - 128.0/3 * eta * (0.5 - xi - eta) * (1 - xi - eta) -  128.0/3 * eta * (0.5 - xi - eta) * (0.75 - xi -eta);
    case 13:
	return 128 * xi * (0.75 - xi - eta) * (1 - xi - eta) - 128 * xi * eta * (1- xi - eta) - 128 * xi * eta * (0.75 - xi - eta);
    case 14:
	return 128 * xi * (xi - 0.25) * (1 - xi - eta) - 128 * xi * eta * (xi - 0.25);
    case 15:
	return 128 * xi * (eta - 0.25) * (1- xi - eta) + 128 * xi * eta * (1 - xi - eta) - 128 * xi * eta * (eta - 0.25);
    }
    return 0;
}

double Triangular_4_Element::det_Jacobi(double xi, double eta)
{
    double J11 = 0,J12 = 0,J21 = 0,J22 = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	J11 = J11 + _dofslist[i-1][0]*phi_xi(xi, eta, i);
	J12 = J12 + _dofslist[i-1][1]*phi_xi(xi, eta, i);
	J21 = J21 + _dofslist[i-1][0]*phi_eta(xi, eta, i);
	J22 = J22 + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    }
    return (J11*J22 - J12*J21);
}

double Triangular_4_Element::xi_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_eta(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_4_Element::xi_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_eta(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_4_Element::eta_x(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][1]*phi_xi(xi, eta, i);
    return -a / det_Jacobi(xi, eta);
}

double Triangular_4_Element::eta_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
	a = a + _dofslist[i-1][0]*phi_xi(xi, eta, i);
    return a / det_Jacobi(xi, eta);
}

double Triangular_4_Element::phi_x(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_x(xi, eta) + phi_eta(xi, eta, i)*eta_x(xi, eta);
}

double Triangular_4_Element::phi_y(double xi, double eta, int i)
{
    return phi_xi(xi, eta, i)*xi_y(xi, eta) + phi_eta(xi, eta, i)*eta_y(xi, eta);
}
std::vector<double> Triangular_4_Element::gradient(double xi,double eta,int i)
{
    std::vector<double> gradient({phi_x(xi,eta,i),phi_y(xi,eta,i)});
    return gradient;
}
double Triangular_4_Element::Global_x(double xi, double eta) 
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][0] * phi(xi, eta, i);
    }
    return a;
}

double Triangular_4_Element::Global_y(double xi, double eta)
{
    double a = 0;
    for (int i = 1; i <= n_Dofs(); i++)
    {
	a += _dofslist[i-1][1] * phi(xi, eta, i);
    }
    return a;
}

int Triangular_4_Element::NdIdx(int i)
{
    
    return _dofslist[i - 1].GlobalIndex();
}

#else

#endif
